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ABSTRACT 
This thesis investigates the use of combined Wavelet decomposition and Wiener 
filtering for the removal of noise from underwater acoustic signals. Several 
Wavelet/Wiener based denoising techniques are presented and their performances 
compared. Performances of the denoising algorithms are compared to those of Wiener 


filter and wavelet thresholding implementation and demonstrate that Wavelet/Wiener 





based methods are also a viable tool for the denoising of acoustic data under more 


restrictive conditions. 
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I. INTRODUCTION 


Underwater acoustic signals are generally embedded in ambient and self noise. 
Ambient noise arises from all natural and man-made acoustic sources and is independent 
of the means used to measure it [1]. It exists in a particular region of the ocean and 
accounts for shipping traffic, marine life, wave motion, wind and numerous other 
sources. The second type of noise, self-noise, is generated by the receiving system, which 
in general consists of a passive sonar receiver and its platform. Usually self-noise can be 
minimized but ambient noise can not be completely avoided. Since noise degrades the 
collected data, it is desired to remove the noise from the contaminated signal before 
further processing, such as detection, estimation or classification, is conducted. 

Detection of signals in noise is a well-researched subject. A Wiener filter is the 
optimum linear filter for detecting the presence of a deterministic signal in noise and can 
be designed to estimate any signal, given that the user has the noise characteristics [2]. 

Wavelet analysis has received increasing interest over the past two decades and 
found wide-spread use in signal processing applications. This is due to its multiresolution 
capabilities that permit analysis of data at different frequencies with multiple levels of 
resolution. Furthermore, wavelets provide the flexibility of choosing a particular wavelet 
function for a specific application. A large variety of orthogonal basis functions can be 
selected for the transformation of a signal. Results have shown that denoising is often 
much easier in the wavelet domain than in the time domain or other traditional transform 


domains. In fact, wavelet noise removal has been shown to be superior to techniques that 








remove noise by conventional filtering. Conventional or convolutional filtering broadens 
edges and hence it makes it difficult to retain sharp features of the original signal [3]. 

The primary goal of this thesis 1s to provide a new approach to the problem of 
removing additive noise from a given signal. In particular, this work will investigate the 
performance of a combined wavelet and Wiener filtering technique to denoise underwater 
acoustic eae: 

This thesis consists of five chapters. Chapter II discusses the theory of wavelet 
analysis by salou with the more familiar Fourier analysis. Chapter III gives a brief 
introduction to Wiener filtering and wavelet denoising. Chapter IV explains the 


algorithms and provides simulation results using synthetic data. Chapter V provides the 


summary, conclusions and proposes some extensions. 














Il. WAVELETS 


Time-frequency signal representations map a one-dimensional time signal into a 
two-dimensional function of time and frequency. Thus, they display time domain and 
frequency domain information. This combination yields a more revealing temporal 
localization picture of the spectral components of a signal. 

The work in this thesis is mainly based on wavelets, which is a very young and 
robust time-scale processing tool. Traditionally, Fourier analysis has been used as a 
robust processing tool in many signal processing applications. Thus, this chapter will first 


discuss Fourier analysis and then present wavelet analysis. 


A. FOURIER ANALYSIS 


Wavelet analysis can be viewed as an extension of Fourier analysis. Thus, this 
section first reviews the main definitions and properties of the Fourier and short-time 


Fourier transforms. 
1. Fourier Series 


Let x(t) be a periodic signal with a fundamental period To (1.e., fo=J/To is the 
fundamental frequency). Then x(t) can be represented as an infinite sum of periodic 


complex exponential functions. The Fourier series representation of x(t) is given by: 


x(t)= Yarn , (2.1) 
a, = os [x@e7 amt dy n=O021L2..:- -. (2.2) 
0 











Equations (2.1) and (2.2) are called as the synthesis and analysis equation, respectively. 
The Fourier series coefficient a, provides the frequency domain information of the signal 
at the n™ harmonic of fo. The coefficient ap is the average value of x(t) over one complete 


period and is given by: 


a, =— [x@ae, (2.3) 


0 T 


which is often referred to as the dc term of x(t). 
2. Fourier Transform 
The Fourier transform of a continuous aperiodic signal is defined by: 
X(f)= [xe at, (2.4) 
which is referred to as the analysis equation. The inverse Fourier transform is given by: 
x)= fx(Ne™*a 25) 


which is the Fourier synthesis equation. Equations (2.4) and (2.5) are known as the 
Fourier transform pair. X(f) gives the decomposition of x(t) in terms of complex 


exponential functions of different frequencies. 
3. Short-Time Fourier Transform 


As can be seen from Equation (2.4), the analysis coefficients X(f) are computed as 
inner products of x(t) with sinusoidal basis functions of infinite duration. The result iS 
that the Fourier transform is a quite robust method for the time-frequency analysis of 


natrow band signals (e.g., sinewaves). It is not matched to non-stationary signals, that is, 











signals whose frequency structures evolve in time. It lacks time localization of 
frequencies, which is especially very important for the analysis of signals of short 
duration. For non-stationary signals, it is necessary to obtain an accurate time-frequency 
representation of the signal that is a function of both time and frequency. By contrast the 
Fourier transform is a function of frequency only. 

The short-time Fourier transform (STFT) was first proposed by Gabor [4]. In this 
transform an analysis window is chosen and this window divides the signal into short 
time segments where the portion of the signals are assumed to be stationary within the 
segments. Initially, the window is located at the beginning of the data allowing the 
computation of the Fourier transform for the first segment. Then, the window is shifted in 
time to process other segments of the data. 


The STFT representation of a signal x(t) is: 
S(t, f)= J x(t)w' (t-te dt, (2.6) 


where w(t) is the analysis window function. The spectrogram of x(t) which provides the 
measure of the signal energy in the time-frequency surface can be obtained from S(7,f) by 
taking its squared magnitude. The time and frequency resolution of the time-frequency 
surface are determined by: 


[e2poceyf at 
At? ==—___—__,, (2.7) 


oo 


Joo dt 


—co 


and by: 











r 2 
[Pwo a 
Af * ==, (2.8) 
2 
[lw cn ar 
where W(f) is the Fourier transform of w(t). 

Equations (2.7) and (2.8) indicate that two impulses in time can be separated only 
if they are more than Af apart, and two tonals must be more than Af apart in frequency in 
order to be discriminated. The drawback of the STFT is related to the Heisenberg 
principle (i.e., uncertainty principle) which is a resolution problem [4] and states that 


there is a natural limit to the resolution, which is given by: 
1 
Ataf =—. (2.9) 
An | 


Using a window of infinite length results in a Fourier transform that provides perfect 
frequency resolution but no time information. Using a short window will give better time 
resolution but poorer frequency resolution. Since the same analysis window is used for all 
frequencies, the time-frequency resolution is fixed (i.e., the duration of all basis functions 
are exactly equal) and the STFT has a time-frequency plane with a uniform tiling [5], as 


shown in Figure 2.1 and Figure 2.2. 


Frequency 


yy Time 


Figure 2.1: Uniform tiling of the time-frequency plane by the STFT with associated basis” 
functions. After Ref. [5]. 
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| ke 
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Figure 2.2: Time-Frequency plane for the STFT; (a) Narrow window, (b) Wide window. 





Due to a constant-bandwidth filter operation, the time-frequency plane is 
partitioned into equal tiles [5]. The analysis filters have identical frequency responses [6]. 


Figure 2.3 shows the frequency response of the STFT. 


Magnitude 


Frequency 


Figure 2.3: Frequency domain of the STFT. From Ref. [6]. 


When the window w(t) is Gaussian, the representation is called the Gabor 
transform [4]. The Gabor transform is optimal in terms of resolution since the Gaussian 
window has the best simultaneous time-frequency resolution for a given window length. 


This meets the lower bound with equality, that is: 





Ataf = m | (2.10) 
Ar 
The discrete Fourier transform (DFT) of a sequence x(n) takes the form: 
N-1 . ee 
X(k)=¥ x(nje ™ k =0,1,2,....N—1, (2.11) 
n=0 : 
and its inverse is as follows: 
aie jzin 
x(n)=—)) X (ke * n=0,1,2,...,N—1. (2.12) 
N t= | 
The discrete STFT is expressed as: 


S(e?”,m) = ¥ x(n)w(n —m)eI™ | (2.13) 


wu=-co 


where m is the center of the analysis window w(n). The frequency variable @ is 


continuous, and takes the range -7<@<S 7. 
B. WAVELET ANALYSIS 


In the last decade, wavelet theory has become an important research area in many 
different fields of mathematics, physics and engineering. Advantages offered by wavelets 
include their efficient coding ability, their localization in time and frequency and their 
computational properties. Some results also claim a somewhat faster performance than 


that obtained by the Fast Fourier Transform. In this section, we will provide an 


introduction to wavelet analysis. 











1. The Continuous Wavelet Transform 


As described earlier, the STFT cannot provide a good resolution in time and 
frequency simultaneously. In other words, abrupt signal changes can only be detected 
within the chosen window length, as the time-frequency resolution is fixed. 

If the signal of interest is composed of high frequency components of short 
duration plus low frequency components of long duration, then the multiresolution 
analysis, which analyses the signal at different frequencies with different resolution, is 
the proper tool. It provides good time but poor frequency resolution at high frequencies 
and good frequency but poor time resolution at low frequencies while the product At Af 
still satisfies the uncertainty principle [6]. 

The continuous wavelet transform (CWT) provides an alternative to the time- 
frequency information of the STFT as it decomposes a signal into a family of two 
parameter basis functions called wavelets. These functions include short duration-high 
frequency and long duration-low frequency components. Thus, time and frequency 
resolution variations are allowed although the time-frequency resolution product, given in 
Equation (2.10), remains constant. Each wavelet basis function is constructed from the 
same function, the so called analyzing wavelet or the mother wavelet. 


The CWT of a signal x(t) with respect to the mother wavelet y(t) is given by: 


1 r * t =f. 
CWT (t, a) =—= | x()w (—dt, 2.14 
re y (— (2.14) 
where T and a are translation and scaling parameters, respectively. The factor I/Va isa 
normalization factor so that the energy of the scaled wavelet is the same as that of the 


mother wavelet. It can be seen that Equation (2.14) is an inner product form, which is 


9 











commonly used to measure the correlation between the data segments and the basis 
functions. 

The mother wavelet y(t) must satisfy certain mathematical conditions. It must 
have compact support, be finite in duration, and have zero mean which requires that it 
must be oscillatory. The name wavelet, which means little wave, comes from its 
oscillation property [1,6]. The scaling factor a is directly related to the frequency 
resolution and inversely related to the time resolution of the CWT. Thus, as a increases 
(i.e., frequency decreases), the wavelet function becomes more stretched in the time 
domain and hence more compressed in the frequency domain. Similarly, as a decreases 
(i.e., frequency increases), y(t) becomes more compressed in the time domain and thus 
more stretched in the frequency domain. This is shown in Figure 2.4. Due to this 
behavior, the CWT can be interpreted as a filter bank composed of band-pass filters with 
proportional bandwidth. 

Thus, the terminology time-scale plane is preferred for the CWT, since the scaling 
factor a in wavelet analysis plays an analogous role to the inverse of the frequency in 
Fourier analysis. 

The squared modulus of the CWT is a distribution of the signal energy in the 
time-scale plane and is called the scalogram. As illustrated in Figure 2.5, the scalogram 
has high time resolution at high frequency (low scale, small a) and high frequency 
resolution at low frequency (high scale, large a). The STFT has a constant time frequency 


tiling as seen in Figure 2.1. 
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Figure 2.4: Scaling and shifting operations of wavelet transform in the time domain and 
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corresponding magnitude behavior in the frequency domain. 
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Time 


Figure 2.5: CWT Time-Frequency plane and associated basis functions. After Ref. [6]. 
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A better understanding of high multiresolution capability of the scalogram 
compared to the spectrogram may be obtained from Figure 2.6. As can be seen from 
Figures 2.6a and 2.6c, the time localization for an impulse function at t=t9 improves for 


low scales in the scalogram while it has a constant region due to the constant length of 





the analysis window in the spectrogram. For a signal composed of two sinusoids of 
frequencies fp and 2fo, the frequency resolution decreases with a in the scalogram as 


shown in Figure 2.6d. 
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Figure 2.6: Comparison of spectrogram and scalogram for two signals. After Ref. [6]. 

















2. The Discrete Wavelet Transform 


The scaling and translation parameters in the CWT can take on any value. Thus, 
this kind of processing is redundant and computationally expensive due to the large 
‘amount of data involved. Furthermore, discretizing is a requirement to perform any 
transform on a computer for real-world applications. This section will provide a brief 
treatment of the discrete wavelet transform (DWT). 

In general, when the energy of the signal is finite and an appropriate wavelet is 
used, not all values of the CWT decomposition are required to reconstruct the original 
signal exactly. Thus, a continuous-time signal can be represented by the knowledge of the 
discrete transform. Although continuous analysis is often easier to interpret, discrete 
analysis ensures space-saving coding and is sufficient for the synthesis [7]. 

The discrete form of the CWT is simply obtained by sampling the time-scale 
plane (i.e., by sampling 7 and a) so that instead of choosing every possible scale and 
position, it uses only a subset of scales and positions at which calculations are made. 


Hence, for a discrete signal x(n) Equation (2.14) takes the form: 


n-b 
a 


), (2.15) 





DWT (a,b) =—= )_ x(n)y *( 
Jaren 

where n, b and a are the discrete versions of ft, t and a, respectively. The scaling factor a 

is discretized first at a=a) ; J=0,1... and the time parameter is discretized with respect 


to the scaling parameter at b=ka= ka, (k is an integer). Thus, the translation steps 


become proportional to the time scaling. This operation yields the decimated DWT at 


scale J and shift k. 
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Setting adg=2 produces scales and positions, which are expressed as powers of 
_ two, and are called dyadic scales and positions. This is the most popular decomposition 
and is used in the standard DWT. This restriction makes the DWT analysis not only much 
more efficient but also provides stable reconstruction, very small reconstruction errors 


and fast algorithms [1]. Hence, Equation (2.15) can be rewritten as: 
DWT (2’ ,k2°)= Er ew’ (2"n-k)=W,,. (2.16) 
2° on | 


An efficient way to implement this scheme using digital filtering techniques 
yields the so-called fast wavelet transform [7] and is known as Mallat's algorithm [8]. It 
uses a highpass filter (HPF) h;(n) and a lowpass filter (LPF) ho(n) which are also called 
quadrature mirror filters (QMF) where: 

h,(n) =(-1)"h,(P-1-n), (2.17) 

H(z) +|H, (2 =1, (2.18) 
where P is the order of the aie sized filter. The filters are complementary filters which 
divide the frequency axis into two equal bands [0-7/2] and [7/2-z] radians, where 7 is 
half sampling frequency, as shown in Figure 2.7a. The signal is decomposed into 
different frequency bands by applying successive highpass and lowpass filtering 
operations. Then the filtered signal is decimated by two to remove the redundant 
information. | 

The output of the LPF branch is called the approximation, and alia the high 


scale, low frequency components of the signal given by: 


yo(k) = > x(n)hy (2k —n). _ (2.19) 











The output of the HPF branch is called the detail of the signal, and contains the 


low scale, high frequency components given by: 


y(k) = ¥ x(n)h,(2k—n). (2.20) 


The decomposition for one level is illustrated in Figure 2.7b. 


LPF approximations 






LPF HPF y(n) 
/\ h,(n) a details 
0 n/2 — y() 
(a) (b) 


Figure 2.7: (a) Spectral partitioning for QMF filters (b) One-level DWT decomposition. 


As can be seen from Figure 2.8, the subband coding BiGuess is repeated on the 
LPF output using the same types of HPF and LPF for further decomposition of the signal 
until two samples are left. Thus, a data of length 2’ is apconpasel into J scales (Mallat 
algorithm). | 

At every level (scale), the width of the lower band is halved by the LPF (i.e., its 
frequency resolution is improved by a factor of two) and its time resolution degrades by 
two, due to decimation by two [6]. Figure 2.9 shows the resulting spectral partitioning. 
The frequency components that are most dominant in the signal will have large DWT 
coefficients, which are located at the scale which includes die Particular frequencies. The 


scales with small DWT.coefficients can be discarded without major loss of information. 
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x(n) > h(n) 4) 


Figure 2.8: Multiple level decomposition. 


Magnitude 


Frequency 


Figure 2.9: Spectral partitioning for the DWT. From Ref. [6]. 


High frequencies will need more samples while fewer samples are required for 
low frequencies. The time resolution decreases as the frequency decreases but the 
frequency resolution increases. In other words, we have good time resolution at high 
frequencies and good frequency resolution at low frequencies. 

The ieee of the original signal x(n) from its DWT coefficients is called 
the synthesis operation. It is performed by following the above decomposition process in 
reverse order as illustrated in Figure 2.10. The signals at every level are first interpolated 
by two (i.e., a zero inserted between each point), passed through the synthesis filters and 


then added. Thus: 
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a(n) = ¥ Ly (ke)hg (2k —n) + y, (kh, (2k -n)]. (2.21) 


k=—co 


The synthesis filter impulse responses are the time-reversed versions of the analysis filter 


impulse responses such as fo(n )=ho(P-1-n) and f}(n)=hi(P-1-n). 


approximation 


coefficients x(n) 


detail 
coefficients 





Figure 2.10: Signal reconstruction. 


In general, the analysis section of the DWT introduces aliasing due to the 
decimation operation. Some amplitude and phase distortions are also possible during the 
analysis step. Therefore, for perfect reconstruction and alias cancellation, the synthesis 
filters must meet certain mathematical conditions and be adapted to the analysis filters 


[9]. 


A two channel filter bank yields perfect reconstruction when [9]: 

F,(z)H,(z)+ F(2)H,(z) = 22", | (2.22) 

F,(z)H,(—2) + F,(2)H,(-z) =0, (2.23) 
where J is the delay in the z-domain. Equation (2.22) is required for no amplitude and 
phase distortion and Equation (2.23) is required for alias cancellation. Thus, choosing 
F,(z)=H,(-z) and F;}(z)=-H)(-z) satisfies Equation (2.23), and removes aliasing [9]. 
Equation (2.23) and a case study related to alias cancellation will be discussed further in 


Chapter IV. 
17 
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Il. SIGNAL DENOISING 


The denoising schemes investigated in this work involve the Wiener filtering of 
the DWT coefficients. The noise removal performance of each scheme is compared to 
that of the traditional Wiener filter and to wavelet thresholding. This chapter briefly 
discusses the Wiener filtering and reviews wavelet-based denoising. It also defines the 


characteristics of the noise model used in this study. 


A. WIENER FILTERING 


In a general linear signal estimation problem, which is illustrated in Figure 3.1, 
the Wiener filter is the optimum linear filter for estimating a signal in additive noise [10]. 
We consider here the finite impulse response (FIR) form of the Wiener filter of length P. 

The noisy observation sequence is represented by: | 

x(n) = s(n) +1(n) : = GB.1) 
where s(n) is the deterministic signal and 7(n) is the noise. 


desired signal 
d(n)= s(n) 










x(n) = s(n) +7(n) 
: Filter 





7 (7) E(n) = s(n) — S(n) 


noise error signal 


Figure 3.1: General linear signal estimation problem. 
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The mean-square error (MSE) between the signal and the estimate for the signal is given 
by: 

o? = E\le(n)l f= E\|s(n)— 5(n)¢ (3.2) 
where S(n)is the output of the filter. The filter criterion selected is the minimization of 


the MSE using the statistical characteristics for the signal and noise [10]. Usually, these 
statistics are not known a priori and must be estimated from the observed signal x(n). 


For a stationary process the output of the filter can be written as: 


S(n) = F ADx(n —l), (3.3) 


1=0 
where h(1) is the impulse response of the filter. Consequently the solution is obtained by 


solving the Wiener-Hopf equation, which has the form: 
P-1 a : 
YR, U-DA (DY =R@ : i=O,L..,.P-1  , (3.4) 
[=0 


where R,(1) is the auto-correlation function of the observations, “*” denotes conjugation, 





and R,,(1) is the cross-correlation function between the desired filter output s(n) and the 
filter input x(n) [2]. If the signal and noise are uncorrelated, assuming a zero mean noise, 


it can be shown that: 
R,O=R,O+R, OD, (3.5) 
R,O=R, 2), (3.6) 


and: 


R,D=R,O-R,O, (3.7) 





where R,(J) and R,(l) are the auto-correlation functions of the signal and noise, 
respectively. An estimate for R,(1) can be obtained from the received data prior to the 


onset of s(n). The MSE is: 
| P-| . 
o2 =R,(0)- YA ORO. (3.8) 
I=0 


It is easier to derive the Wiener-Hopf equation using matrix notations. In this case, 


the filter output is expressed as: 
P-i 
$(n) =" hOx(n—) =[X)Vh, ) = 469) 
l=0 


where h is defined as: 


h(0) 
h(1) 
he| = f. (3.10) 
h(P-1) 
and: 
x(n—P +1) x(n) 
x(n—P) x(n) 
x(n) = : x(n) = (3.11) 
x(n) x(n—-P +1) 


The matrix form of the Wiener-Hopf equation takes the form [2]: 
Rh’ =F... (3.12) 
In the stationary case, R* =R, and hence Equation (3.12) can be written as: 


R h=r (3.13) 


x SX * 


The solution for the optimum filter coefficients is: 
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x “sx? 


h,,, =Ry fF, 
while the MSE is given by: 

o2? =R.(0)-F7h’ =R,(0)—-hF_, (3.15) 
where R, = E4x(n)x"™ (n)} and r= E{s(n)ix(n)I }- When the signal and noisé are 
uncorrelated and have zero mean then R,=R,+Ry. Thus, the signal correlation matrix R, 


can be determined from R,=R,-R, and rf, can be extracted from the first column of R,. 


When the signal s(n) is non-stationary, the Wiener filter becomes time-varying 
and take on the form h(n,1) [2]. Although the same equations can be used to determine the 
optimum filter weights, a short-time Wiener filtering P gorithm is a more efficient 
approach to this problem [11]. 

The Wiener filter procedure is illustrated in the following graphs. Figure 3.2 plots 
the application of the Wiener filter to a noisy sinusoid of length N=1024 with a constant 
frequency of f=0.05 at an SNR level of -3 dB. A Wiener filter length of P=30 was 
chosen. The top three pists display the original, noisy and denoised signals, respectively. 
The Wiener filter, as shown in the bottom plot of the figure, enhances the frequency 
bands where the signal level is strong and attenuates the frequency bands where noise is 
present. Note that the frequency domain plots have a normalized frequency axis. 

The performance of the Wiener filter in terms of the normalized MSE is displayed 
in Figure 3.3 for various combinations of filter and data lengths, normalized signal 
frequency, and SNR levels. The normalized MSE that will be used for all comparisons is 


defined as: 


E= ; Jae (62) : 3.16 
ceria —_ 





























where N, s(n) and $(n) respectively represent the data length, the noise free and denoised 


signal. Thus, the MSE chosen compares the signals which are energy normalized. 
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Figure 3.2: Denoising with the Wiener filtering. Sampling frequency f;=1]. Left column: 
time domain signal. Right column: spectrum in dB. 
The top plots display MSE versus SNR levels between —15 to 15 dB. The bottom plots 
display MSE versus the normalized frequency spanning the range from 0 to 0.5. The 
graphs represent the average MSE values obtained from ten independent trials. From the 
left-hand side plots of this figure, it is clear that the MSE decreases as the filter length 


increases. In addition, the performances of the filter orders 12, 20 and 30 are close to 


each other. 
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A data size of 1000 or more points was used to obtain reliable estimates of the 
correlation matrices [1]. However, the right-hand side plots of Figure 3.3 indicate that 
increasing the length of the observed sequence does not significantly improve the MSE. 
As a result, it seems reasonable to choose P=/2 and N=1024 from the viewpoint of the 
computational cost factor. Thus, these values will be used later in the performance 


comparisons of the denoising algorithms investigated in this study. | 


f=0.3 , N=1024 | f=0.3 , P=12 





| 
| 
SNR (dB) . SNR (dB) 
SNR=0 dB , N=1024 SNR=0 dB , P=12 





(0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 
~ Normalized Frequency Normalized Frequency 


Figure 3.3: The performance of the Wiener filter. The left two plots are for various filter 
lengths. The right two plots are for various data lengths. 








B. WAVELET DENOISING 


As mentioned earlier, wavelet signal decompositions have numerous applications 
in signal processing including noise suppression. It has been proven that denoising in the 
wavelet domain is often much easier than in the other conventional domains. Moreover, it 
preserves the sharp features of the signal to be recovered, unlike the traditional 
techniques which use convolutional filtering [3]. 

Wavelet shrinkage, introduced by Donoho and Johnstone [12,13,14], is a 
nonlinear processing technique that can be employed for data compression or denoising 
according to the thresholding method used. Recall that in the wavelet decomposition, 
each transform coefficient can be viewed as a measure of the correlation between the 
signal and decomposing basis. Thus, the transform will contain a few large coefficients 
for good correlation and many small coefficients for poor correlation. Although all the 
coefficients are required for perfect reconstruction, the general shape of the original 


signal can be represented by a small number of coefficients carrying significant signal 


_ energy (i.e., coefficients with large magnitudes). 


Wavelet shrinkage algorithm keeps only those important coefficients based on a 
nonlinear thresholding. It is a discarding operation of the coefficients that fall below a 
given magnitude. Wavelet denoising consists of three steps [7], namely: 

1- Taking the wavelet transform of the data by decomposing it into J levels with 

a suitable wavelet basis function, 
2- Selecting a threshold T for each level eit 1 to J and soft or hard thresholding 


the wavelet (detail) coefficients, 


pa 





3- Performing inverse wavelet transform to reconstruct an estimate of the 
underlying signal using the nema approximation coefficients and the 
modified detail coefficients. 

Hard thresholding zeroes out the coefficients W),; of magnitudes lower than the 


threshold value and is given by: 


wi = fe Wal >7 (3.17) 


0, We | ST 


Soft thresholding is an extension of the hard thresholding is given by: 


WS, = nats XW, a -T), W, | >t (3.18) 


0, W,| ST 
It is clear from Equation (3.18) that soft thresholding first zeroes out the coefficients of 
magnitudes smaller than the threshold value, and next shrink the magnitude of the 
remaining coefficients towards zero. Note that W,’, has discontinuities while W;,, does 


not. 


Both thresholding methods mentioned above may be used in denoising and data 


compression. In general, hard thresholding is used in data compression while soft 


thresholding is usually employed for denoising. 

The noisy data can be expressed as: 

x(n) = s(n) +ow(n), (3.19) 
where s(n) is the signal to be recovered from the noisy data, w/( n) is zero mean additive 
white Gaussian noise (AWGN) and o is the standard deviation of the AWGN. Since the 
orthogonal wavelet transforms are linear operations, the wavelet coefficients will be in 


the form of Equation (3.19). 


26 





In the J-level wavelet decomposition of x(n), the approximation coefficients at 


every successive level become less and less noisy. The high frequency components in the 
signal due to the noise show up in the detail coefficients and mostly at the finest scale [7]. 
The finest scale is the frequency band between 7/2 and z, where 7 corresponds to the half 
sampling frequency. The goal of the wavelet-based denoising is to retain the coefficients, 
which contain signal information while removing those due to the noise contribution. 
Thresholding is a lossy operation and usually provides a lowpass version of the original 
signal. Thus, it cannot provide perfect reconstruction of the original signal since 
thresholding the detail coefficients will also remove some of the signal components [9]. 
Algorithms computing the threshold value T require an estimation of o. Several 
variants include Stein’s Unbiased Risk Estimate (SURE), fixed threshold, Minimax 
criteria or a combination of those. The Matlab Wavelet Toolbox from MathWorks used in 
this study employs four thresholding methods, in its routine function thselect, as follows 


[7]: 


1- Rigrsure: Adaptive threshold selection using SURE, 

2- Sgtwolog: Fixed _ threshold equal to 2log(N) ; 

3- Heursure: A mixture of the previous options, 

4- Minimaxi: Fixed threshold for achieving minimax performance in terms of the 

MSE, 

where N is the number of wavelet coefficients. 

Results show that wavelet denoising performances improve as the SNR level 
increases and the frequency decreases. We also noted that thresholding schemes fail at 


SNR levels 0 dB and below. 


27 











The MSE performance comparisons of the threshold selection methods of the 
toolbox are displayed in Figure 3.4. The test signal is a cosine wave of 1024 points 
duration with a constant frequency of 0.3 at a sampling frequency of 1. As mentioned in 
the DWT section of Chapter II, this signal can be decomposed into 10 levels. However, a 
5-level wavelet decomposition covers almost the complete frequency spectrum except 
very low frequencies. The sym8& (i.e., Symlet8) wavelet was empirically chosen and the 
signal was decomposed into 5-level wavelet coefficients. For every level, a threshold 
value calculated according to the given threshold selection method and soft thresholding 
was applied to the coefficients. The resulting MSE indexed by SNR is shown in the top 
part of Figure 3.4. The bottom part of the figure plots MSE versus the normalized 
frequency at a fixed SNR level of 0 dB. All points represent averages of 10 independent 
runs. Without loss of generality, it can be seen that the accuracy of all four methods 
decreased as the frequency increased. As expected, the error of all techniques decreased 
as the SNR increased. From both plots, it is clear that the fixed form threshold has the 
worst performance. The reason for this might be due to the fact that the threshold 
obtained by sgtwolog was set too large to keep all the signal information coefficients. On 
the other hand, the rigrsure and heursure methods attained the best performance for all 
conditions. 

The soft and hard thresholding method: are compared in Figure 3.5 where the 
rigrsure wavelet thresholding method is applied to the same noisy tonal used in the 
previous section. The signal was decomposed into 5 levels using the sym8 wavelet. We 


note that soft thresholding leads to a better performance in the MSE sense. 
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Figure 3.4: Comparison of different threshold selection methods. 


We noted that wavelet denoising schemes usually threshold the detail coefficients 
only, as done for example in the Wavelet toolbox. However, we also noted that, n cases 
where the signal ‘ryosiaation is not contained in the approximations, denoising 
performances were improved by thresholding the approximation coefficients, too. This is 
to be expected as in such cases the approximation coefficients contain noise only, and as 


a result should be disregarded. 
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Finally, some knowledge about the signal in advance is useful in noise reduction 
using wavelet thresholding, since choosing a suitable wavelet basis function that matches 
the signal character will provide a few but large magnitude wavelet coefficients that 


represent good correlation. 
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Figure 3.5: Denoising with the wavelet thresholding. Sampling frequency f,=J. Left 
column: time domain signal. Right column: spectrum in dB. 
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C. NOISE DESCRIPTION 


It was desirable to make a few assumptions of the additive noise model used to 
provide simplicity in the study and analysis of the denoising algorithms that we 
investigated in this thesis. 

Results have shown that the ocean acoustic noise is usually a colored Gaussian 
random process [1]. However, in this study it was assumed to be zero mean additive 
white Gaussian noise (AWGN) based on the fact that some pre-whitening techniques can 
be employed if necessary. This assumption was very helpful to allow a wide variety of 
simulations for different underwater signals ‘al comparisons using different filtering 
techniques. | 

The noise was assumed to be a stationary process. This assumption holds as 
generally denoising algorithms are concerned with short durations of data of a few 
seconds. In such short periods of time, the properties of the ocean itself and the ambient 
noise change very little [1]. 

Next, we assumed that the variance of the AWGN was unity and varied the signal 
amplitude to obtain various SNR levels. This also enabled us to make comparisons of our 
algorithms to those of Donoho’s study. 

Finally, we note that our results are based on the knowledge of the exact noise- 
only sequence. In practical applications, one does not have access to the exact noise 
information which then must be estimated. These estimates are expected to lead to 


performance degradations in the various schemes studied. 
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IV. WAVELET/WIENER-BASED DENOISING 


The previous chapter described how Wiener filtering and wavelet thresholding 
can be applied to remove noise from acoustic signals. In this chapter, we will present four 


denoising algorithms which combine the DWT and Wiener filter. Our goal is to 





investigate whether these combinations can achieve a better performance than those 
obtained with either the classic Wiener filter or wavelet thresholding in a MSE sense. The 
main idea behind the schemes investigated is to apply FIR Wiener filtering to the data 


transformed to the wavelet domain. 


A. DENOISING USING ORTHOGONAL DECIMATED WAVELET 


TRANSFORMS 


Consider the following 2-channel multirate system, which is essentially a QMF 
analysis/synthesis bank with additional Wiener filtering operations, represented by Wo 


and W, inserted in each bank. Recall that the original QMF setup is alias-free. 


n(n), x(n) 





Analysis Decimators Wiener Expanders Synthesis 
Filters Filters Filters 


Figure 4.1: Decimated Wavelet/Wiener-based denoising, /-level decomposition. 
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The system input is the noisy signal x(n icidadla ), where s(n) and 7(n) are the 
deterministic signal and the noise, respectively. Given that we have access to a noise only 
segment of noisy data, the one-level DWT is first applied to this noise-only segment. 
Next, the DWT of the noisy data is computed. The correlation function of the wavelet 
filtered noise and noisy signal determine the characteristics of the Wiener filter. The 
outputs of the Wiener filters represent the modified approximation and detail scehticients 
which will be used for the reconstruction to get ii estimation of the clean signal s(n), as 


shown in Figure 4.2. 


n(n) 


x(n) s(n) 





Analysis Decimators Wiener Expanders Synthesis 
Filters Filters Filters 


Figure 4.2: Wavelet/Wiener-based denoising scheme. 


The denoising operation can be easily expanded to multiple levels using the 
orthogonal DWT decomposition. Figure 4.3 illustrates the J-level Wavelet/Wiener 


filtering procedure. 
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n(n), x(n) ® 7 Wz) 
@a}e z 
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Figure 4.3: The analysis step of the J-level Wavelet/Wiener-based denoising. 


Note that in a perfect reconstruction QMF bank, the reconstructed signal is merely 
a scaled and delayed version of the input. However, introducing Wiener filters in the 
above system produces aliasing and distortion, as the analysis/synthesis relationship is 
modified. Note that aliasing is undesirable and must be canceled to achieve a good 
estimation of the signal s(n). The next sections will present different approaches in order 
to overcome this problem. 

Figure 4.4 depicts the application of the Wavelet/Wiener scheme to the noisy 
sinusoid already considered in the previous chapter. A Wiener filter length of 12 and a 5- 
level decomposition is used. The first two top plots show the noisy signal and the 
resulting Wiener filtered signal. Denoised signals obtained with the combined method are 
shown in the third and fourth plots, where dbJ0 (..e., Daubechies10) and db40 wavelets 


are used, respectively. 
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Note that choosing a higher order wavelet reduces aliasing effects and the MSE, 
as well. Moreover, the Wavelet/Wiener denoising gives a better performance than the 
usual Wiener filtering in a MSE sense. However, they also both show the presence of 


aliasing in the denoised signal. 
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Figure 4.4: Denoising with the Wavelet/Wiener method. Sampling frequency is 1. Left 
column: time domain signal. Right column: spectrum in dB. 
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Figure 4.5 displays the spectrum of the tonal at normalized frequency equal to 0.3. 
For this frequency the Wiener filter gives a MSE of 0.1833. Although both spectra 


indicate aliasing, they have still better MSEs than those obtained using the Wiener filter. 


MSE=0.1226 


Magnitude (dB) 
Magnitude (dB) 





o.1 oO.2 o.3 u 
Normalized Frequenc 


Figure 4.5: Denoising with the Wavelet/Wiener method using db10 and db40 wavelets. 5- 
level wavelet decomposition. 


B. DENOISING USING 1-LEVEL ORTHOGONAL UNDECIMATED 


WAVELET TRANSFORMS 


We showed in the previous section that introducing the frequency band Wiener 
filters resulted in creating aliasing in the denoised signal. This section investigates a /- 
level high-pass/low-pass (i.e., band-specific) Wiener filter denoising scheme in the 
wavelet domain, where the decimation and interpolation operations have been removed, 
as illustrated in Figure 4.6. The main idea is to compare the band-specific Wiener filter to 
the classic Wiener filter in the absence of further distortions created by the 
decimation/interpolation operations present in the decimated wavelet transform. 
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Simulations showed that this implementation introduces amplitude scaling in the 
filtered signal due to successive filtering operations. To allow visual comparison we 


normalized the output signals to have unit energy. 


n(n), x(n) s(n) 






Figure 4.6: Undecimated Wavelet/Wiener-based denoising. 


Figure 4.7 illustrates the noisy, band-specific Wiener filtered signals using two 
different wavelets (i.e., db10 and db40) and the respective spectra obtained for the single 
tone at normalized frequency 0.05 used previously in Section A. The length of the 
‘Wiener filter is 12 and SNR=-3 dB. The plots show that no aliasing or distortion is 
present. As can be seen the MSE values are higher than those of the scheme with 
decimators and expanders, however, they are still lower than that of the classic Wiener 
filter. 

Figure 4.8 compares the spectra of the denoised signals obtained with the Wiener 
filter, the decimated Wavelet/W iener, and the undecimated Wavelet/Wiener denoising 
schemes using a single sian at normalized frequency equal to 0.3 and an SNR level of —3 
dB. The Wiener filter length is 12 and the signal is decomposed into 5 levels in the 
decimated Wavelet/Wiener filter scheme. 

Results show that increasing the wavelet order in the wundecimated 


Wavelet/Wiener scheme does not seem to improve the MSE performance Significantly. 
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Figure 4.7: Undecimated Wavelet/Wiener Denoising scheme. Left column: time domain 
signal normalized to unit energy. Right column: spectrum in dB. 
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Figure 4.8: Filtering scheme comparisons. One single tone at normalized frequency f=0.3 
and SNR=-3 dB. 
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C. DENOISING USING ORTHOGONAL DECIMATED WAVELET 


TRANSFORMS WITH ALIAS CANCELLATION 


We showed in previous sections that aliasing present in the decimated — 
Wavelet/Wiener filter scheme is a direct result of the decimation/expansion operations 
present in the decimated wavelet transform. In this section, we consider two schemes 
designed to cancel the aliasing present in the denoised signal when applying the 
decimated transform. The first one attempts to modify each wavelet synthesis filter 


independently of the others, while the second one uses cross-band filters. 


1. Alias Cancellation Synthesis Filters 


Consider the following J-level decimated Wavelet/Wiener denoising scheme 


where the filters G;(z) (k=0,/) are introduced for alias cancellation. 





n(n), x(n) 





Figure 4.9: J-level decimated Wavelet/Wiener filter scheme with alias cancellation filters. 
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The filters G;(z) are to be defined so that the denoised signal 5(m) contains no 


aliasing term. 

The z-transform of the outputs to the analysis filters Hp and H; are given by: 

To(z)=Ho(z)X(z); 

Tj(z)=Hi(z)X(2).- (4.1) 
The z-transforms of the decimated signals u;(n) are given by: 

Uo(z)=1/2[Holz'7)X(2"” + Hol -2' )X(-2 1, 

U;(2)=1/2[H(27)X(27) +H )X(-2")]. (4.2) 
The z-transforms of the Wiener filters' outputs v;,(n) are: 

Volz)=1/2[ Wolz)Holz"7)X(z"”)+ Wolz)Hol-2'”)X(-z" II, 

Vi(2)=1/2[W(2)H (27 Xe” )+ Wil ZH -2)X(-2")). (4.3) 
After the interpolation steps, we have: 

Yo(z)=Vo(2)=1/2[ Woz )Ho(z)X(z)-+ Wo(z Hol -z)X(-z)], 

¥j(2)=Vi(2)=1/2[Wi( 27 )H(2)X(z)+ Wiz JH (-2)X(-2)]. (4.4) 
Next, the z-transforms of the terms z;(n) are given by: 

Zo(2)=1/21Fo(z)Wolz’ )Ho(z)X(2)+Fo(z) Wo( 2’ JHol-z)X(-2)], 

Z(2=1/2LF (2) Wi C)H DX) +F (2) Wi -2)X(-2)). «SS 


Finally, the z-transform of the reconstructed signal 5(7) is: 

S$ (z) =Golz)Zol(z)+G(z)Z1(2). (4.6) 
Substituting Equations (4.1) to (4.5) into Equation (4.6) leads to: 

§(z) =1/2X(z)[Golz)Fo(z) Wo( 2 )Ho(z)+ Gi(2)F (2) Wiz )H(2)) 


+1/2X(-z)[ Go(z)F oz) Wo”) Hol -z) + GZ) F (z)Wi 2 )H(-z)). (4.7) 
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Recall that canceling the alias term X(-z) is required to get perfect reconstruction. 
It can be seen from Equation (4.7) that the aliasing term can be canceled by choosing the 
filters such that the term Go(z)Fo(z)Wo(2)Ho(-z)+Gil2JF (z)Wi(Z)H1(-z) is equal to zero. 
Therefore, aliasing is canceled by selecting the filters Fj(z) and G;(z) as: 

Fo(zj=Hi(-2)  , = F'i(z)=-Hof-2), (4.8) 

GolZJ=Wiz) , — Grlz)=Wol2’). (49) , 


The QMF bank perfect reconstruction property already implements Equation (4.8). 





- Equation (4.9) can be implemented after designing the Wiener filters. 


Replacing Fo(z), F(z), Go(z) and G;(z) by their values in Equation (4.7) leads to 
S(z) =1/2X(2)[Wi(2)H1(-z)Wol 2 )Ho(z)-Wolz7)Ho(-z)Wi(2)H,(z)]. Using the fact that 


H,(z)=-2%Ho(-z') leads to §(z) =Wo(2)Wi(z)z™X(z). 





Simulations results showed that this scheme did not give satisfactory results as it 
introduced further distortions in the denoised signal, even though the aliasing terms tend 


to be canceled 15]. Thus, the second alternative considered is presented next. 
2. Alias Cancellation Analysis Filters 


In this section, we consider a decimated Wavelet/Wiener denoising scheme, 
_ Shown in Figure 4.10, where K;(z) (k=0,/) are alias cancellation filters whose lengths are 
the same as those of the Wiener filters W,(z). As done earlier, the alias cancellation filters 
will be derived by using the z-transform. 
The z-transform of the outputs to the analysis filters Hp and H are given by: 
To(z)=Ho(z)X(Z), 


T}(z)=H)(z)X(z). | (4.10) 
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Figure 4.10: J-level decimated Wavelet/Wiener filter scheme with alias cancellation 
analysis filters. 
The z-transforms of the decimated signals u,(n) are given by: 
Uo(z)=1/2[ Holz” )X(z'")+ Hol-2'”)X(-2"" 1, 
U(2)= 1/22” )X(2"7) + Hy -2"7)X(-2')]. (4.11) 
The z-transforms of the Wiener filters outputs v;{n) are: 
Vol(z)=1/2[ Wo z)Holz’”)X(z'?)+ Wo(z)Ho(-z")X(-2"), 
Vi(z)=1/2[Wi(z)Hi(2'7)X( 27 )+ Wiz) (2 )X(-2" Nh. (4.12) 
The z-transforms of the alias cancellation filters outputs z;(n) are: 
Zo(z)=1/21 Kol 2) H(z )X(z'7) + Kol 2) H(-2'”)X(-2' 
Z1(z)=1/21K (2H 2” JX (27) + Kz) Hol 2" )X(-2'")]. (4.13) 
After summation and interpolation steps, we have: | 
Yo(z)=Vol2)+Zolz = 1/2. Wol 2 )Holz)X(z)+ Wolz )Hol-z)X(-2)] 
$1/2[ Ko 2 )H(z)X(z)-+Ko( 2 A -2)X(-2)], 
Y(2)=Vi(2)+Z(Z )=1/21Wi(z’)H (z)X(2z)+W(z )i(-2)X(-z)] 
+1/2[K (2 )Ho(z)X(2z)+ KZ )Hol-2)X(-2)). (4.14) 
Finally, the z-transform of the reconstructed signal S(7) is: 
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S(z) =Folz)¥o(2)+F (2)¥ (2). (4.15) 
Substituting Equations (4.10) to (4.14) into Equation (4.15) leads to: 
S(z) =1/2X(z)[Fo(2)Wo(< Hol z)+F (z)WiZ H(z) 
+Fo(2)Ko(Z Ai(2)+F (2)K iz Hol 2)1, 
+1/2X(-2)[Fo(z)Wolz’ )Hol-2)+F 1(2)Wi(z )Hi(-2) 
+Fo(z)Ko( )Hi(-2)+Fi(z)Ki(¢ )Ho(-2)]. (4.16) 
It can be seen from Equation (4.16) that the aliasing term X(-z) can be canceled 
when: 
Fo(z)Wo(z’ )Ho(-z)+F 1(z)Wi(Z’ )Hi(-z)+ Foz) Kol 2 )Hi(-2)+F (z)Ki(2’ Hol -z)=0. 
(4.17) 
Recall that Fo(z)=Hi(-z) and F}(z)=-Ho(-z) can be selected to insure the analysis and 
synthesis filters properties lead to perfect QMF reconstruction. Substituting these 
properties into Equation (4.17) and rearranging leads to: 
Fo(2)F (z)[Wilz")-Wolz 1+ Fo(z) Kol2’)-F (2) Ki(z )=0, (4.18) 
and, 
Fo(z)Fi(z)[Wilz)-Wolz N=Filz)Ki(<)-Fo(z) KolZ). (4.19) 
For a wavelet filter of order of L (i.e., filter length is L+1), Wiener filters and alias 
cancellation filters of order of P (ie., filter length is P+1), the filters 


Fo(2),F (z),Wo(Z),Wi(Z),Ko(Z) and K}(z’) are given by: 


Fo(z)=agtajz'+...4azZ™. (4.20) 
Fi(z)abotbiz'+..4b12. | (4.21) 
Wo(z2)=mo+m)z?-+...+mpz??. a (4.22) 
Wiz )=notnjz’+...+npe?”. (4.23) 

















Kol 2 )=yotyiZ +... tYpZ (4.24) 
Ki(2 Jatotz1e + AZPe (4.25) 
Thus, 
Fo(z)Fi(z)=aobot (agb;+ajbo)z" +...rarbiz 
=hothjz'+...+hgyz, (4.26) 
where hp=aobo, hj=aob)+ajbo,..., haz=axDx. 
F(z)°=Fo(2)Fo(z)=(agtajz'+...+412"). (aptaj)z' +. tazz") 
=Apdot(apaj+ajag )z? +... apart =kotkiz! +... re oe aaa (4.27) 
where kp=aodo, kj=A0a)+a1,..., K2p=azaz. | 
FA(2P =F (2)F (2)=(botbiz'+...+b1z"). (bot bjz'+...+biZ") 
=bobo+(bob +b bo) +... tb z= got B12 t--.+8k”, (4.28) 


where go=bobo, g j=bob ] +b bo,.. -s S2L=OD 1. 


Thus, the left hand side (LHS) of Equation (4.19) may be rewritten as: 


Fo(2)F AZ) Walz )-Wolz )1=Uo(no-mo)] + Uhi(no-mo)]z'+...+[har(np-me)]e™. 


(4.29) 
Next, using matrix notations, the coefficients for each z factor in the LHS of Equation 


(4.19) may be written as: 


h oO. 0 h, O 0 
h 0 : h 0 
> Ay > hy 
hy: | ; Ng hy, My 
0 > td 2 jt 0 a ed! Ss Sy (4.30) 
hy, 0 Np > hy 0 Mp 
0 hy | (+1) : 0 hy | (+11 
O: as Thy O°. QO? ve. Mey 
[2¢L+P) +1 {P+1] 


{20L+P)+1 >{P+1] 
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which may be rewritten as: 


h OO . O}h O 
h O > |h, 0 
: Oh, ile. . Me 
hy, > yy 

OO: “a. | oe 

hy, Ol: hy 

0 hy}: 0 

0 0 h,|0 0 
[2(L+P)+1 >{2(P+1)] 





(4.31) 








{2(P+1)x1] 


Recall that the right hand side (RHS) of Equation (4.19) is expressed as: 


Fi(2)Ki(2)-Fo(z)’ Koz’). 


(4.32) 


Expanding F';(z)’Kj(2’) and Fo(z)’Ko(Z) leads to: 


F(z)’ Ki(2)=go0z0+(g120)Z' +...+(g21zp)z 
Fo(z)’ Kol 2 )=koyot(kryo)z'+...+(koryp)z2 


Thus, the RHS of Equation (4.19) becomes: 


eee (4.33) 


(4.34) 


F(z) Ki(2)-Fo(z)’Ko(Z)=(goz0-koyo) +(g 1Zo-kyo)Z' +...+(gorzp-koryp)z 


(4.35) 


In matrix notation, the coefficients for each z factor in Equation (4.35) can be written as 


g oO .. O 
§ 0 
. 80 
Sor: Xo 
0 : 
$21 0 ip 
0 2, | (eH) 
0 0 $21 
[2(L+P)+1pqP+1] 


0 0 
O . 
K | 
. Yo 
a eon oe (4.36) 
Kor 0 yp 
0 ky} (@+xt] 
0 Keay 
(2(L+P)+1Pq¢P+1] 
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and, 
Zo 0 O1k, O 0 
8 0 : | K, : 
&o Ky : 
8a, Ky, | 
0 =O. -< ae = (4.37) 
0 
Sa 0 Koy 0 . 
0 0 
; - 
° ° ° ; ° [2(P-+1)x1] 
0 0 a. ay 0 0. se. Koy ‘ 
(2(L+P)+1)¢2(P+1)] 
A 


Thus, Equation (4.19) can be expressed in matrix notation as: 

Ax=B. (4.38) 
Note that the matrix equation to be solved has a larger number of mine than columns. 
Therefore, the problem can only be solved in a least square sense. Thus, solving Equation 
(4.38) indicates that the aliasing will be minimized but not exactly canceled. Therefore, 
the least square solution is obtained by: 

x=B\A, | (4.39) 
where the vector x provides the coefficients for the filters Ko( z) and K,(z). 

Figure 4.11 illustrates the spectra of the noisy and the denoised signals obtained 
with the Wiener filter, the decimated Wavelet/W tener Ge the scheme of section A 
without alias cancellation filters) and the decimated Wavelet/Wiener with alias 


cancellation analysis filters denoising schemes. The test signal is a sine wave with a 


length of 1024 samples at normalized frequency 0.2 and SNR level —3 dB. The length of 


the Wiener filter is 12 and the wavelet function is sym8. For the decimated 


Wavelet/Wiener denoising scheme of section A, a J-level decomposition is used, as well. 
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The plots show that aliases are canceled using the scheme with alias cancellation filters. 


As can be seen that the MSE value is also lower than that of the classic Wiener filter. 


Figure 4.12 displays a portion of the time domain signals and Figure 4.13 shows 


the frequency responses of the band-specific Wiener and alias cancellation filters of the 


denoising scheme. 
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Figure 4.11: Algorithms comparison; One tone at frequency f=0.1, sym8. 
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Figure 4.12: Time domain signals. 
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Figure 4.13: Wiener and alias cancellation filter characteristics. 
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Figure 4.14 shows that when using a db/0 wavelet function on the same signal as 
considered in Figure 4.11, the aliased ‘ain located at normalized frequency f=0.3 is not 
canceled. Our simulations showed that the amount of aliasing cancellation depends on the 
specific selection and combination of Wavelet and Wiener filter used. 


Expansion of the scheme to multiple levels is left as a further study. 
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Figure 4.14: Algorithms comparison; One tone at frequency f=0.1, db10. 
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D. DENOISING USING NON-ORTHOGONAL UNDECIMATED WAVELET 


TRANSFORMS 


This section considers the application of the non-orthogonal wavelet 
decomposition to the denoising of underwater signals. Recall that sampling the CWT in 
powers of two leads to the standard DWT (1.e., dyadic scales and positions) which is 
orthogonal, decimated and time-variant due to the decimation. Although those linear 
transforms are computationally efficient, design of analysis and synthesis filters is 
constrained by the orthogonality restriction, that is, the filter impulse responses are 
directly related to each other. However, non-orthogonal transforms have fewer | 
restrictions on the decomposition filters, thus, better frequency resolution thus can be 
achieved than with orthogonal transforms. 

The analyzing wavelet function used in the non-orthogonal transform is 
implemented using the undecimated A-trous al pai and a modified version of the 
original Morlet wavelet [16,17] and given by: 

9 ae A 


Bet 


wth=e"e 7? , (4.40) 
where 77 represents the center frequency of the mother wavelet and B represents the roll- 
off parameter, which is a constant proportional to the bandwidth of the analyzing 
wavelet. The Fourier transform of Equation (4.40) takes the following form: 
1 (o-1)" 
w(@)=—V2me 7 . (4.41) 
p 
It is clear from Equation (4.40) that one can change the 7) and f parameters of the 


resulting decomposition, which results in different spectral partitioning. 
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The frequency resolution of the non-orthogonal transform may be modified by 
introducing subfilters (called voices) defined in each scale. These voices allow the user to 
vary the spectral partitioning of the transform [17]. Note that the partitioning of the 
orthogonal decompositions is fixed when a basis function is chosen. Thus, using non- 
orthogonal decomposition, much better narrowband frequency partitioning can be 
achieved than with orthogonal decompositions. However, a possible drawback with 
voices is the loss of temporal resolution due to the uncertainty principle, and an increase 
of computational cost per octave proportional to the number of voices selected. More 
details on analysis and synthesis operations of the non-orthogonal wavelet transform can 
be found in [17,18]. 

In this scheme, using the analyzing wavelet shown in Equation (4.40) will lead to 
complex DWT coefficients. Moreover, although those eeeticienls are associated with the 
white noise input, they are correlated unlike those of orthogonal transforms, which are 
uncorrelated. | 

Figure 4.15 display the spectral partitioning obtained by the A-trous 
decomposition with various roll-off parameters and number of voices for a 3-scale 
decomposition using a center frequency of 0.857. Figure 4.16 depicts the 3-level spectral 
partitioning with different number of voices and center frequencies for a fixed roll-off 
factor of 0.15. In order to have a full coverage of the frequency axis, it is necessary to use 
a high roll-off parameter, which leads to larger overlap between the subfilters, as shown 
in Figure 4.15. Moreover, to cover very low frequencies one should use a high 
decomposition level. Finally, note that changing the center frequency 7 of the mother 


wavelet shifts the sub-filters in the spectrum, as shown in Figure 4.16. 


52 














2 voices — beta: 0.15 


Magnitude 
» 8 Oo ® O 
0 0 0 0 90 0 
) 

P 


2 voices — beta: 0.35 


0 


Magnitude 
N @Q 


= 
K 


0 
ity 


0.2 0.4 0.6 
Normalized Frequency 





4 voices — beta: 0.15 













; 
Ho | 






Figure 4.15: Spectral partitioning obtained for the A-trous algorithm, 7=0.857. 
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Figure 4.16: Spectral partitioning obtained for the A-trous algorithm, B=0. 15. 











The denoising performance of the non-orthogonal undecimated transform is 
illustrated in Figure 4.17. We applied a 5-scale decomposition to the noisy sinusoid, 
which is at normalized frequency 0.05 as shown saree The signal consists of 1024 
samples and SNR level is —3 dB. We used 20 sub-filters in the spectrum (i.e., 4 voices per 
octave), a roll off factor 0.25 and a center frequency 0.757. The resulting MSE is 0.1032. 
Recall] that the classic Wiener filter gives a MSE of 0.1994 for a filter length of 12. 

Simulations showed that this implementation introduces amplitude distortion like 
the scheme investigated in section B. Thus, we normalized the signals to unit energy for 


ease of comparisons. 
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Figure 4.17: Denoising with non-orthogonal undecimated wavelet transform. 


54 





























E. METHOD COMPARISONS 


In this section, two synthetic signals, Doppler and Decay, are studied to compare 
the different methods. These signals are of interest to us since they have some essential 
features of real world acoustic signals. The signals, consist of 1024 data points and have 


an SNR of —3 dB, as shown in Figure 4.18 and Figure 4.19, respectively. 
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Figure 4.18: The clean and noisy Doppler, SNR=-3 dB. 
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Figure 4.19: The clean and noisy Decay, SNR=-3 dB. 
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The denoising schemes compared and parameters used for Doppler are: 


1- 


Oe 


Wiener filter with a length of 12, 

Wavelet thresholding with 5-level sym8& decomposition using soft heursure 
thresholding method, 

aliaied Wavelet/Wiener filtering with 5-level sym8 decomposition using a 
multiresolution Wiener filter length of 12, 

Cudecinated Wavelet/Wiener filtering with J-level sym8& decomposition 
using a band-specific Wiener filter length of 12, 

Decimated Wavelet/Wiener filtering with alias cancellation filters (i.e., K(z) 
where k=0,1) and J-level db10 (since it performed better than sym8 for the 
Doppler test signal) decomposition using a band-specific Wiener filter length 
of 12, 

Non-orthogonal undecimated wavelet transform with 5-level decomposition 
and a Wiener filter length of 12, 2 voices per octave using a roll-off factor B 


of 0.35 and a center frequency 77 of 0.857. 


The resulting denoised signals and MSEs obtained for the energy normalized 


signals are displayed with the order shown above in Figure 4.20. Note that the 


undecimated scheme leads to amplitude distortion in the denoised signal as expected. On 


the other hand, the A-trous algorithm does not have an amplitude distortion but a delay 


unlike the situation in the previous section. Moreover, it has a worse MSE value than any 


other schemes due to the delay. As a general comment, it can be said that all denoising 


schemes lead to similar results except the A-trous algorithm. In addition, it is clear from 


Figure 4.20 that there is difference between visual and numerical quality of denoising. 
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Although the A-trous, thresholding and decimated Wavelet/W ones denoising schemes 
achieve the highest MSEs, they provide a smoother visual reproduction of the clean 
signal. 

Next, we applied the denoising schemes mentioned earlier to the second test 
signal Decay. For the wavelet thresholding scheme, different wavelet filters were tried. 
Since db10O gave better results, it was used instead of sym8. In addition, a roll-off factor 
of 0.25 and a center frequency of 0.957 were used for the non-orthogonal undecimated 
wavelet denoising. All remaining parameters are the same as those of the Doppler signal 
experiment. 

Figure 4.2] depicts the denoised signals and the corresponding MSEs. As 
expected, the undecimated schemes have mapieade distortion. In addition, wavelet 
thresholding significantly performed worse than the other schemes. 

Note that knowledge of the exact noise-only segment was assumed in the 
simulations conducted above as mentioned in Chapter Ill, section C. Results showed that 
performance degradations occurred in the denoising schemes in absence of this 
assumption, as expected [15]. 

Finally, we note that one can make performance improvements by choosing 


parameters based on experimental evidence. 
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Figure 4.20: Performance comparisons of the denoising algorithms. Note the disparity 
between MSE and visual quality. 
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Figure 4.21: Performance comparisons of the denoising algorithms. 
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V. CONCLUSIONS 


In this thesis, we investigated a number of denoising schemes by combining 
wavelets and Wiener filtering techniques. The performance comparison of those schemes 
to that of the classic Wiener filter and wavelet thresholding was conducted by applying 
them to the problem of acoustic noise removal. 

In general, simulations show that the Wiener filter with an appropriate length 
outperforms other denoising techniques in terms of MSE, especially at lower SNR levels. 

The combined Wavelet and Wiener filter denoising scheme also performs well 
under more restrictive conditions. However, the performance is highly dependent on 
signal of interest, accurate estimation of noise and wavelet basis function used. Aliasing 
is another problem with this type of schemes due to the introduction of Wiener filters, 
although they provide low MSEs. 

Two types of alias cancellation filters were conducted: single-band post synthesis 
and cross-band filters. Results showed that the single-band post filters reduce the aliasing 
component but introduce some additional distortion in the denoised signal. Cross filter 
transfer functions were obtained as the result of solving a least square matrix equation 
which showed that this option leads to minimization of the aliased components but not to 
their exact cancellation. Results showed this scheme performance to be quite sensitive to 
the choice of wavelet filter type, Wiener filter order and signal frequency location. Thus, 
considering the added problem generated by the alias canceling filters, users should 
consider using the denoising schemes without any added alias canceling filters when 


some amount of aliasing is permissible. 
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Finally, the undecimated A-trous Wavelet/Wiener denoising scheme was 
considered. Simulations showed better performance results than those obtained with the 


classical Wiener filter, however, the computational cost is much higher. 
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